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Abstract 

Epileptic seizures are one of the most well-known dysfunctions of the nervous sys- 
tem. During a seizure, a highly synchronized behavior of neural activity is observed 
that can cause symptoms ranging from mild sensual malfunctions to the complete loss 
of body control. In this paper, we aim to contribute towards a better understanding 
of the dynamical systems phenomena that cause seizures. Based on data analysis and 
modelling, seizure dynamics can be identified to possess multiple spatial scales and 
on each spatial scale also multiple time scales. At each scale, we reach several novel 
insights. On the smallest spatial scale we consider single model neurons and investigate 
early-warning signs of spiking. This introduces the theory of critical transitions to ex- 
citable systems. For clusters of neurons (or neuronal regions) we use patient data and 
find oscillatory behavior and new scaling laws near the seizure onset. These scalings 
lead to substantiate the conjecture obtained from mean-field models that a Hopf bifur- 
cation could be involved near seizure onset. On the largest spatial scale we introduce 
a measure based on phase-locking intervals and wavelets into seizure modelling. It 
is used to resolve synchronization between different regions in the brain and identifies 
time-shifted scaling laws at different wavelet scales. We also compare our wavelet-based 
multiscale approach with maximum linear cross-correlation and mean-phase coherence 
measures. 

Keywords: epileptic seizure, multiple time scales, multiple spatial scales, wavelets, crit- 
ical transitions, excitable systems, phase- locking, correlation measures. 
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1 Introduction 



Trying to predict epileptic seizures using time series analysis has been an important research 
topic for decades. In particular, the now wide-spread use of EEG (electroencephalography) 
techniques to acquire data has been a major driving force of the subject. The review article 
[55] and the recent book [76] provide perspectives what has been achieved in seizure pre- 
diction. The main goal was to identify and characterize a pre-ictal phase occurring before 
the onset and to design measures that approximately predict the critical starting time of 
the seizure [19]. Since research has focused in this direction there are still gaps [69] in our 
understanding of seizures from a dynamical systems perspective [861 US] . In this paper, we 
are going to address this issue and focus on dynamical mechanisms [TS] instead of aiming at 
a predictive technique for seizures. 

The main themes of our results are the deep links to mathematical multiscale techniques 
[68| [65] and the observation of scaling laws at different spatio-temporal levels. From models 
based on biophysical principles of brain dynamics it is expected that multiple spatial [10] and 
multiple time scales [32] play an important role for epileptic seizures [67] . Based on a combi- 
nation of analyzing epileptic seizure patient data and neuron modelling we split the problem 
into three spatial scales and show that at each individual spatial level the problem exhibits 
multiple time scale behaviour. We point out that our approach to verify the existence of 
multiscale phenomena is primarily data-driven and complements modelling approachs (see 

e.g. m- 

On the smallest spatial scale, we employ model-based analysis of single neurons [51} [57] 
using a multiple time scale stochastic FitzHugh-Nagumo model [221 EOl HZ] with a focus 
on early-warning signs [75j, scaling laws and control failure of spiking. In particular, we 
investigate three different cases of spiking and provide the first results of critical transition 
theory [39] for neurons in an excitable state. Critical transition theory for systems without 
equlibria near bifurcations has recently been applied successfully in climate modeling ^441 
[1] and in ecological systems [151 HI] but apparently has not been applied to neuroscience 
problems yet. We analyze three different regimes for the relationship between noise and time 
scale separation and show that the variance can be a precursor of spiking in some parameter 
regimes while it fails in the low noise case. In this context, we point out that the distributions 
of interspike intervals [IB] has been studied extensively in single neuron models but that our 
work only studies the time series locally near a bifurcation and does not require multiple 
events. 

The second spatial scale which we consider are clusters/regions of neurons [631 [77]. Here 
we use electrocorticogram (ECoG) data; see Appendix [X] We examine the onset of the 
epileptic seizure using the variance as a simple univariate measure. We observe that during 
a certain period before the seizure the variance shows oscillations. Furthermore, very close 
to the transition to a seizure the inverse of the variance displays a linear scaling law. Based 
on critical transition theory, these observations are generically characteristics for Hopf bifur- 
cation [10] . It is very important to note that many seizure models [72l [T!] [501 ESI [9] suggest 
a Hopf bifurcation as a main mechanism as the transition point. Therefore, our results not 
only provide a first application of local scaling laws near bifurcations to data but also vali- 
date the proposed bifurcation mechanism arising from biophysical principals. Similar to the 
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individual neurons scale, we point out that distributions of interseizure intervals have been 
studied [H3] but that we do not require multiple events. 

On the largest spatial scale we analyze the synchronization and correlation between differ- 
ent brain regions [43j . Several bivariate measures have been proposed [55j to study epileptic 
seizures but the underlying complex network structure makes the problem difficult [H] . Our 
approach utilizes a recent technique calculating phase-locking intervals (PLIs) [38] based on 
wavelet transforms [HS] . Wavelet-based methods have been applied previously in the context 
of epileptic seizures [S] but our approach is the first to investigate PLIs and associated phase- 
locking. We show that our wavelet-based method [381 ES] measures increasing phase-locking 
and resolves a multiple time scale structure near the seizure onset. Furthermore, we observe 
a linear scaling law of average phase-locking and that phase-locking at different scales often 
starts at different times where the lower scales tends to increase earlier compared to higher 
scales. These results apply near the seizure onset and could potentially relate to recently 
observed rapid discharges [871 ES] • We also compare our results to other bivariate measures 
such as maximum linear cross-correlation [731 l2l] and mean phase coherence [H]. 

In summary, our study introduces two recently developed methods (critical transitions, 
wavelets/PLIs) into the analysis of epileptic seizures. Using critical transitions theory we 
give the first analysis of early-warning signs for excitable neurons, identify a potential Hopf 
bifurcation as the seizure onset mechanism from data and find a new scaling law of single- 
event time series data at the cluster level. For the wavelet-based phase-locking technique, we 
provide a comparative study to other bivariate measures and discover a scaling law occurring 
at time-shifted onset times. On each of the three spatial levels we were also able to identify 
a multiple time scales structure, based on a data-driven time series approach. 

2 Single Neurons 

We start on the level of single neurons. Clearly it is very problematic to get data in this 
case before epileptic seizures so that we resort to model neurons. The main question will 
be whether we can predict a spike in the voltage time trace of the model neuron before it 
occurs. The FitzHugh-Nagumo (FHN) model [22] EOl EO] is a simplification of the Hodgkin- 
Huxley equations [31J which model the action potential in a neuron. We point out that the 
methods we are going to present here are going to apply to a much wider class of excitable 
neuronal models than the FHN equation such as the original Ho dgkin- Huxley model [71] 
or the Morris-Lecar system [27] since these models have similar bifurcation structure and 
multiple time scale properties [31]. 

There are several forms of the FHN-equation [28]. One possible version suggested by 
FitzHugh is the Van der Pol-type [18] model 



where x represents voltage, y is the recovery variable and 7, 6, e are parameters. We think 
of b as an external signal or applied current [48j and assume that the time scale separation 



dx 



ex = X — x^ — y, 
y = ^x-y + b, 




(1) 
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e satisfies < e ^ 1 so that x is the fast variable and y the slow variable. The dynamics 
of ([T]) can be understood using a fast-slow decomposition [TSl [25] • Setting e = in ([T]) 
gives a differential equation on the slow time scale r defined on the algebraic constraint 

Co:={{x,y)eR'^ ■.y = x- x^ =: Co(x)}. 

We call Co the critical manifold; see Figure [U Differentiating y = co{x) implicitly with 
respect to r we find y = x{l — 3x^) so that the differential equation on Cq can be written as 

■yx — X + x^ + b 



which we refer to as slow fiow. Observe that the slow fiow is not well-defined at the two 
points {x±,y±) = {±x^^^^ , co{±x^^^^)) . Applying a time re-scaling to the fast time t := r/e 
to (lH) gives 




Figure 1: Simulation of ([3]) with e = 0.005 and 7 = 2 using an Euler-Maruyama numerical 
SDK solver pQ]; red curves are deterministic trajectories with a = and blue curves are 
sample paths with a = 0.0028. Systems have always been started at {xo,yo) = (0,-0.2). 
The critical manifold Co is shown in grey and the ?/-nullcline as a dashed black curve, (a) 
6 = 0, the equilibrium for the full system lies on Cq. (b) b = 0.8, the equilibrium lies on 
Cq~ near the fold point The deterministic trajectory has only one spike while 

noise-induced escapes produce repeated spiking for the stochastic system. 

Setting e = in [2] gives the fast fiow where y' = implies that y is viewed as a parameter 
in this context. Observe that Cq consists of equilibrium points for the fast fiow and that the 
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points {x±,y±) are fold (or saddle-node) bifurcation points [79] in this context. The critical 
manifold naturally splits into three parts 

Cq~ := Co n {x < Cq := Co H {x_ < x < x+}, Cq^ := Cq n {x > x+} 

where Cg^ are attracting equilibria and Cq are repelling equilibria for the fast flow. We view 
Cq~ as the refractory state and Co"*" as the excited state for the neuron. For e = trajectories 
are concatenations of the fast and slow flows. We will consider two different situations for 
the parameters (7, b). In the first situation we chose the parameters so that ([1]) has a single 
equilibrium point on Cq fl F where T := {{x,y) E M."^ : y = '~fx + b} is the y-nullcline of the 
FHN-equation; see Figure [T]^al)-(a2). For e = suppose that {xo,yo) := {x{0),y{0)) G Cq~; 
then the slow flow moves the system to a jump via the fast subsystem to C^'^ 

occurs, the slow flow on Cq^ brings the system to {x+,y+) and another jump returns it to 
Cq^ . This is the classical relaxation oscillation [23 However, in neuroscience one often 
also considers the excitable regime [M] where the global equilibrium {x*,y*) for the system 
is stable and lies on Cq~ close to ?/_); see Figure [T]^bl)-(b2). In this case, a trajectory of 
([1]) can generate, depending on {xo,yo), at most one excursion/spike to the excitable state 
before returning to {x*,y*). Repeated spiking in the excitable regime can be obtained using 
the more general stochastic FHN-equation 

6Xj- X-j- Xj. y ~\~ (^^Ti /Q\ 

Vt = 'JXr-yr + b, 

where is delta-correlated white noise (^ri^T2) = ^(^i — ^2) and a is a parameter represent- 
ing the noise level. We can now ask whether individual neuron spiking activity already has 
precursors. This viewpoint should provide new insights how neurons are able to control syn- 
chronization and how control failure occurs. Recent results on predicting critical transitions 
[75] suggest that statistical precursors can be used to predict events similar to spiking in 
neurons from a time series without knowing their exact location. The detailed mathematical 
theory can be found in [39], HQ] . 

Here we present the first application of this theory in the context of single neurons. We 
want to predict a spiking transition from a neighborhood of Cq^ to Cq^ and consider the 
variance as an early-warning sign 

Vg := Var(xs) restricted to Xs near Cq~ . 

Observe that we can view Vs also as a function of y, and write V = V{y), since the 
mapping between yg and s is bijective when restricting to Cq~ . In the relaxation oscillation 
regime (see Figure [T|^al)-(a2)) and if (e, a) are sufficiently small it can be shown \^ S\ that 

V{y) ~ , as 2/ ^ (4) 

Vy - y^,- 

for some constant A > and where \ye,- — y- \ is small. Therefore an increase in fast voltage- 
variable variance can potentially be used to predict and to control spiking if no equilibrium 
exists near (x_, ?/_). Here we extend the results of [40] by investigating the excitable regime. 
Figure El shows an average of the variance V computed over 100 sample paths using a sliding 
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Figure 2: Average of the variance V{y) = Var(a;(y)) (black curves) over 100 sample paths 
starting for s = at {xo,yo) = (—1,0) up to a final time T. The green curves are fits 
of V using (jlj) with fitting parameters A and y^^^. Fixed parameter values are (e,7) = 
(2,0.005). (a) Relaxation oscillation regime with {b,a,T) = (0,0.02^6,0.28). (b) Excitable 
regime with {b,a,T) = (0.8, 0.02A/e, 0.7); sample paths can exhibit oscillations around the 
stable focus equilibrium {x*,y*) which are visible in the variance, (c) Excitable regime with 
{b,a,T) = (0,0.05^6,0.8) where larger noise regularizes the variance similar to (a), (d) 
Excitable regime (6, a, T) = (0, 0.005A/e, 0.9) where smaller noise does not allow fast escapes 
from {x*,y*) and yields decreasing variance. 

window technique [39]. Figure |2]^a) shows the relaxation oscillation regime where we can 
confirm the theoretical prediction (jlj). 

The excitable regime is much more interesting since the equilibrium point {x*,y*) can 
lead to a variety of distinct regimes depending on the noise level. In Figure EJ^b) the noise is 
at an intermediate level so that deterministic oscillations around the equilibrium are visible 
in the variance before an escape; hence the prediction (jlj) is not a good prediction of a spike 
but one should rely on the oscillatory mechanism before escapes. In Figure [2](c) the noise is 
larger which provides a regularizing effect for the variance via noise-induced escapes. This 
relates to the well-known mechanism of coherence resonance [17]. In Figure [2](d) the noise is 
very small so that sample paths need exponentially long times to escape and are metastable 
near {x*,y*). This causes a decrease in variance and will make predictions very difficult. 
The different scaling regimes for noise level and time scale separation are discussed in more 
detail in [201 EH EH HQ] . 

Based on our results we can conclude that control of spiking could depend crucially on 
noise level and statistical properties of the state of a neuron. In particular, in the excitable 
state already a small change in the noise level or system parameters can result in a substantial 
loss of control due to unpredictable spiking. This could cause undesirable synchronization 
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and continuous spiking. Let us point out that this is just one possible explanation for 
a potential prediction/control failure during epileptic seizures but our results show that 
prediction at neuronal level can already be extremely complicated. Therefore we proceed to 
look at the next scale in our analysis and move from single neurons to clusters/regions of 
neurons. 

3 Local Data & Clusters 

On the level of regions, we can start to analyze data obtained before epileptic seizures. The 
eight time series we use are described in detail in Appendix \M A natural extension of 
our previous strategy is to compute the variance for each time series using a sliding window 
technique and to understand the scaling laws associated with the variance on the cluster level. 

Figure [3] shows the results of this computation. We plot the inverse of the variance Vj~^ 
for i G {1, 2, . . . , 8} since this makes it easier to understand the scaling of Vi near the seizure 
point at t = tc- All four plots have some features in common: 

(A) Vj^^ decreases near the seizure. The scaling law seems to be given by 

Vi ~ -j——:, as t tc. (5) 

(B) There are multiple local maxima and minima for V~^ before approaching the seizure 
point. This indicates that we should expect oscillations in statistical indicators near 
epileptic seizures. 

(C) The last local maximum before shows that there is also a period of low variance 
close to a seizure. 

(D) The last local maximum before tc is already very close to the seizure. This means that 
predictions could be very difficult just based on a calculation of the variance. 

The next problem to consider is what types of dynamical models can reproduce the 
behavior we have observed from the data analysis i.e. we look for a model for the variance 
in clusters/regions of neurons that displays the observed oscillatory behavior and scaling 
law. At first glance, Figures [3]^a)-(h) could be interpreted as a summation of Figure [2](b) 
i.e. of neurons that are (almost) in synchrony where the coherent spiking originates from 
the noise-induced escape of a spiral sink. However, the real problem in understanding the 
dynamical mechanism of epileptic seizures is shown in Figure H] where we also plot the inverse 
of the variance V~^ near a critical transition. The similarities to the data in Figure |3] are 
clear; all four observations (A)-(D) also apply in Figure |H The data in Figure H] have been 
generated using a simple model for a Hopf critical transition [23 SO] : 

eil.r = yTXl,r-X2,r + Xi,r{xl^ + xl^^) + a{au^i^^ + ai2^2,T), 

eX2,r = Xi^r + yTX2,T+X2,ri^lr+^2,T)+^i(^21^1,r + a22^2,T), (6) 
Vt = 1, 
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Figure 3: The eight plots show the average channel activity Mj (top, blue) and the average 
of the inverse variance V^^ = 1/Vi (bottom, black) for the eight time series i G {1, 2, . . . , 8}; 
the horizontal axis is the time axis where the labels correspond to the sample point number. 
The green dots mark local maxima of which correspond to local minima of Vi. The 
fitted red curves are linear and demonstrate that the variance increases near the epileptic 
seizure. The black dashed vertical lines are inserted for orientation purposes and roughly 
mark a seizure onset beyond which a very high variance is observed. 



where ^j^r are independent white noise processes that satisfy (0,ti0,t-2) = ^(''"1 ^ '^2) for 
j = 1,2. The model ([6]) was first analyzed in the context of delayed Hopf bifurcation 
[6T| |62]. Observe that the deterministic part of the fast variables {xi,X2) is the normal form 
of a (subcritical) Hopf bifurcation ^42j. The slow variable y can also be viewed as time since 
= [t — To) + I/O- For the simulation in Figure H] we have chosen 

e = 0.0005, a = O.OOlv^, an = 1 = 022, 012 = 0.2 = 021 (7) 

with a deterministic initial condition {xi^Q,X2fi,yo) = (0,0,-0.3). It is known that near a 
sub- or supercritical Hopf bifurcation a scaling law of the form holds ^40j . Obviously the 
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Figure 4: Time series of the fast variable X2 (top, blue) and the associated inverse of the 
variance = Var(x2(y))~^ (bottom, black) for a Hopf critical transition model with 
parameter values [71 cf . also Figure O 

scales differ between Figure [3] and Figure H] but those can be re-scaled to match. Therefore 
we have found a dynamical model that could potentially explain the qualitative features of 
a single variance time series for a cluster of neurons. 

It is very important to note that we have obtained the conjecture that a Hopf bifurca- 
tion is involved in the transition to a seizure without a detailed biophysical model. In fact, 
several mean-field models for various types of epileptic seizures do exhibit Hopf bifurcations 

[HI E21 EH EDj that form a boundary between a regular equilibrium (non-seizure) an os- 
cillatory (seizure) regime. However, there are several mean-field models available [17] and 
also other bifurcation mechanisms have been identified to play a role near seizure onset |84j . 

Our methods also have another important implication regarding the distinction between 
a preictal and a proictal state |81j. From a dynamical perspective, it was suggested that one 
can differentiate between models that show a distinct preictal state with a parameter driving 
the system to a bifurcation or systems showing a proictal state where noise-induced escapes 
play a dominant role [16]. A subcritical Hopf bifurcation is a model that can interpolate 
between the two cases. Consider in the following two cases: 

(1) e = 0, 5 > and y < 0: the equilibrium (xi,X2) = (0,0) is a stable focus for the 
deterministic dynamics but it is well-known [23] that a finite-time noise-induced escape 
always occurs. This can be viewed as the transition beyond a basin boundary given by 
the unstable limit cycles [8^. If we include another (seizure-state) attractor beyond 
this basin boundary we can view the situation near (0:1,0:2) = (0,0) as a "purely 
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proictal" state. It is well-known how to calculate the probabilistic likelihood of this 
Hopf transition and also for many other bifurcations involving met ast ability [21] • 



(2) e > 0, < 5 ^ 1: If the noise is sufficiently small then we will reach the Hopf 
bifurcation point with high probability [5] and our prediction method via scaling of 
the variance and critical transitions applies. We are in a "purely preictal" situation. 

Obviously there is a continuum of possibilities in between these two situations [59l [5] 
depending on the scaling of noise and time scale separation. In fact, the results shown in 
Figure [2] illustrate the variation in such a continuum situation for the saddle- node bifurcation. 
A similar study of intermediate regimes for all bifurcations, including the Hopf bifurcations 
on a mean-field level, can certainly be carried out similar to the strategy employed in |10] . 
This study is beyond the scope of this paper as we focus on the multiscale features and basic 
dynamical analysis techniques here. 

4 Correlations between clusters 

In the preceding two sections we investigated neuronal dynamics at different spatial scales, 
from single model neurons to neurophysiological data from clusters of neurons, using the 
variance as a univariate measure. In the following section, we will focus on the dynamics from 
many clusters of neurons encompassing a larger spatial scale. In contrast to the previous, 
bivariate measures for the activity between different clusters will be used. 

Bivariate measures can take into account the correlation of two signals. Information about 
the correlation of neuronal activity between different anatomical regions can give insights 
into the state of the network as a whole. With regard to epilepsy, correlation based measures 
such as mean phase coherence (MFC) and maximum linear cross correlation (MLCC) have 
yielded promising results in identifying pre-ictal states |571 [SH [771 [21] . 

In this section, we will start by considering the maximum linear cross correlation for the 
ECoG data used in the preceding parts, reviewing and confirming some recent observations. 
We will then continue to extend the bivariate analysis to wavelet-based synchronization 
measures able to resolve pairwise correlations at different frequency bands. We will compare 
these results to those obtained using MLCC and MFC. Our focus is again on multiscale 
character of the system with the goal of identifying scaling relationships at each level of 
observation. 

4.1 Maximum linear cross-correlation 

The maximum linear cross-correlation (MLCC) quantifies the similarity between two time 
series Fi{t) and Fj{t). MLCC is a linear measure of lag-synchronization which captures the 
normalized product of two time series dependent on a lag r [7B] : 



= max 





max 



T 



10 



where 



N-T 



(9) 



t=i 



is the hnear cross-correlation function. As a measure of synchronization between activity in 
different anatomical areas, MLCC has been proposed and successfully applied as a precursor 
for pre-ictal brain activity jMj 156] . We computed the MLCC of 5 randomly chosen signal 
pairs for each time window (5000 sampling steps, a consecutive time window being shifted 
50 sampling steps forward) . Figure [5] shows the average over the 5 pairs for each of the 8 
time series considered in the preceding sections. 
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Figure 5: Maximum linear cross-correlation MLCCi for eight pre-ictal time series i G 
{1,2,...,8}. Vertical lines indicate the approximate onset of the seizure attack. 



In most of the depicted time-series (patients 1, 2, 3, 4, 6, 7) an increase in the MLCC can 
be observed with the seizure onset (FigJS]) which is in agreement with the general observa- 
tion that seizure attacks are characterized by an increased synchronization between cortical 
regions [l3]. Prior to epileptic seizures a decrease in MLCC values has been reported and 
used to identify a preseizure state [Ml ES]. To relate to these reports and later also com- 
pare MLCC to the wavelet-based synchronization measure {PLI) (see following section), we 
calculated MLCC for a pre-ictal and an inter-ictal time interval. Figure M (left column) de- 
picts the time series of MLCC values of patient 4 during a pre-ictal (top) and an exemplary 
inter-ictal interval (middle), an interval being at least 6 hours apart from the next seizure 
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attack. Average values of MLCC are plotted left in the bottom row illustrating the com- 
parably lower values during pre-ictal intervals. MLCC levels are lower during the pre-ictal 
compared to the inter-ictal interval confirming recent reports of decreased synchronization 
as one characteristic precursor for a seizure. 





Figure 6: Decrease of synchronization measures during a pre-ictal interval. Left column: 
time series of maximum linear cross correlation during a pre-ictal (top) and an inter-ictal 
(middle) interval. Right column: time series of (PLI) for three scales during a pre-ictal (top) 
and an inter-ictal (middle) period are depicted. Vertical dashed lines indicate the onset of 
the seizure attack. Averages over the first 150000 sample points of each time series indicate 
a distinct decrease of each synchronization measure during the pre-ictal interval (bottom 
row). 

A lot of effort has been put forward to utilize the observed synchronization drop in 
predicting seizure attacks, most of these works addressing the question whether it could be 
used to identify a preseizure state [53 EH ESI EH [771 EI]. this work we are not addressing 
this issue but focus on the dynamics and scaling relations of correlation measures near the 
seizure onset. For this purpose we extend the analysis to wavelets able to resolve correlations 
between clusters for different frequency bands. 
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4.2 Wavelets 



Wavelet analysis has been applied in neuroscience research for some time [I3l [12]. Wavelet 
coefficients Wk provide a frequency-dependent moving average over a time series which can 
be used to derive a time-resolved frequency-profile for the data given. This capacity has also 
been made use of in the detection of seizures [80] and the invest igationen of frequency profiles 
of epileptic seizures in humans and animals [Bl ES]- Wavelet analysis [65j can also be used 
as an elegant tool to identify intervals of phase synchronization (or phase- locking) between 
neurophysiological time series. The phase definition can thereby be used for broad-band 
synchronization analysis or analysis of a specific frequency of interest. 

In this study, we investigated broad-band phase-locking between pairs of signals as intro- 
duced in ^88j. There, the original signal is decomposed with respect to multiple scales related 
to frequency bands of decreasing size. To derive a scale-dependent estimate of the phase dif- 
ference between two time series, we follow the approach described in [38] using Hilbert 
transform derived pairs of wavelet coefficients [88]. The instantaneous complex phase vector 
for two signals Fi and Fj is defined as: 

_ W,{F.^W,iF,) 

where denotes the k-th scale of a Hilbert wavelet transform and ^ its complex con- 
jugate. A local mean phase difference in the frequency interval defined by the k-th wavelet 
scale is then given by 



A0„(t) = Arg{Q,,), (11) 

with 

{WkiF,yW,iF,)) 



C it) = \-- V >^\-3J/ 

V{\Wk{Fmm{F,)\^) 

being a less noisy estimate of Cjj averaged over a brief period of time At = 2'^8 [38j. 
One can then identify intervals of phase-locking (PLI) as periods when |A0j j(t)| is smaller 
than some arbitrary threshold which we set to 7r/4 here. We denote phase-locking intervals 
between two signals Fi and Fj as PLIi j. To obtain a measure of frequency-specific phase- 
locking in a defined time window, we calculate the sum of PLIi j for all pairs of signals and 
normalize this expression to confine the measure to the interval [0, 1]: 



n 



steps 



where Usignais is the number of signals and nsteps the number of time steps in the time window 
under consideration. 



We analyzed data for each patient for 3 different scales, referring to frequency bands 25- 
12, 12-6 and 6-3 Hz for patients 1-3, 5-8 and 32-16, 16-8 and 8-4 Hz for patient 4, respectively. 
The computation of {PLI) was done for time windows of 5000 sampling steps, consecutive 
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Figure 7: Phase-locking measure {PLI)i for the eight time series i G {1,2,..., 8}. Colors 
correspond to different scales. The vertical dashed lines indicate the approximate onset of 
the epileptic seizure attack. 

time windows were shifted forward by 50 sampling steps. Figure [7| shows the results of this 
computation. 

In all 8 patients, comparably low values of (PLI) {{PLI) < 0.5) are observed for all 
scales. Similarly to MLCC synchronization as measured by the phase-locking intervals for 
different scales is decreased during an pre-ictal interval compared to an inter-ictal one (Fig. [6l 
right column). The observed decrease in synchronization measure suggests that application 
of (PLI) could also prove useful in preseizure state detection algorithms, similar to the 
MLCC. 

As mentioned earlier, our focus is on the dynamical behavior near the seizure onset. Aside 
from the aforementioned low values of (PLI), some characteristic features can be observed: 

(A) Phase-locking measured by (PLI) increases around seizure onset times. (Similar to 
MLCC, this is seen less clearly in patients 5 and 8.) 

(B) The increase of (PLI) for different scales often starts at different times. Phase-locking 
of lower scales tends to increase earlier compared to higher scales. 

(C) The increase of (PLI) appears to be linear. 

Point A comes not as a surprise reflecting the fact of increased synchronization between 
cortical regions observed during seizures. For better clarity regarding point B and C, we plot 
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the inverse of (PLI), (PLI)^^, for times close to the seizure onset (Fig. [8]). In many of the 
patients (1, 2, 3, 6, 7) one can observe that {PLI)~^ of lower scales decreases earlier than 
higher scales meaning that (PLI) of higher frequencies shows an earlier increase compared to 
lower frequencies. In fact, the temporal order in which (PLI) increases often appears to be 
directly correlated to the numerical order of scales (see patient 1 in fig. [8] for example). We 
furthermore observe (PLI)^^ to decrease with j suggesting a linear increase of its inverse 
(PLI) with time. Fitting a linear function (PLI) oc t close to seizure onset times also 
provided the better fit compared to power-law or exponential relationships (Fig. E]). 

Another nonlinear measure based on phase synchronization is the mean phase coherence 
|571 [H] . For two pairs of neurophysiological time series Fi and Fj it is given by 



with A(j)ij{t) being the phase difference between the two signals at time t and () denoting 
the average over time. We calculated Rij for all pairs of signals using the wavelet-derived, 
scale-dependent phase differences for each patient. The average (R) over all Rij showed a 
similar time course as (PLI) (Fig. [9]). Near seizure onset, the same temporal order of the 
increase in synchronization was observed indicating independence from the specific measure 
of phase-synchronization. Direct comparison of both nonlinear synchronization measures (R) 
and (PLI) to MLCC suggests that the frequency resolved measures add new information 
at the onset of the seizure. Therefore such multiscale measures may potentially be better 
suited to explain the dynamical process that causes a seizure attack. 

5 Discussion 

In the present paper we aimed for a better understanding of the dynamical processes in- 
volved in seizure generation. Our approach extended over three spatial scales involving two 
recently developed methods (critical transitions and wavelet derived phase-lock intervals). 
We showed for the first time that the theory of critical transitions HP] can be applied 
in the context of excitable neurons operating near the spiking threshold. On the level of 
clusters of neurons we identified a potential Hopf bifurcation as the seizure onset mechanism 
from data based on this theory and found a new scaling law of single-event time series data. 
On the largest spatial scale we observed a scaling law occurring at time-shifted onset times 
and compared our wavelet-based phase-locking measure to other bivariate measures. 

One of our main results is the observation of scaling laws on different spatial scales - for 
individual neurons @, for activity of clusters of neurons ([5]) and for the increase of phase- 
locking near the seizure onset. A recent publication highlighted five power- law scaling laws 
related to epileptic seizures and their analogy to earthquakes (the Gutenberg-Richter distri- 
bution of event sizes, the distribution of interevent intervals, the Omori and inverse Omori 
laws and the conditional waiting time until next event) |6l]. Other works investigating 
scaling laws of ictal and interictal epochs reported similar inter-seizure-interval statistics in 
genetically altered rats while in human data no power-law distribution was observed [831 IH2] • 




(14) 
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Figure 8: The inverse (PLI)^^ = 1/ {PLI) for the eight time series is plotted near the seizure 
onset. The fitted go with 1/t and demonstrate that the variance increases hnearly near the 
epileptic seizure (black dashed vertical lines). 
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Figure 9: Comparison between (PLI) and mean phase coherence (R) for patient 1. Both 
measures based on phase-synchronization show a similar behavior. Colored vertical lines 
indicate the beginning of the increase in synchronization near the seizure onset (black dashed 
vertical line). Both measures show an earlier increase for lower scales. 



The observation of such scaling laws is important because it may guide new models of 
seizure dynamics by allowing insights into the dynamical processes that may have generated 
the underlying data. Many of the scaling laws reported here and elsewhere [6l] exhibit power 
laws. The observation of similar scaling laws on different spatial scales, from single neurons 
to the size distribution of different seizures, strongly emphasizes the multi-level character of 
epileptic seizure generation. More importantly, it yields insights into the dynamical proper- 
ties of the underlying system [35]. The strong analogies between seismic shocks and brain 
seizures have previously been pointed out and hypothesized to emerge from the structural 
commonality of the two systems: both are composed of interacting nonlinear threshold os- 
cillators and are far from equilibrium [3^ . Critical dynamics is believed to be a consequence 
of these structural properties in both these systems. Recent findings in preparations of 
rat cortex p] and primate brain in vivo [66] exhibiting power-law statistics of activity, a 
hallmark of phase transitions [21 |45l |51] , have led to the hypothesis that also human brain 
dynamics is poised at a phase transition [U [3H]. Although such statistics can result from 
different processes, the self-similar behavior captured by the diverse scaling laws on different 
levels might potentially be related to the notion of criticality in brain dynamics. Models 
describing epilepsy should also resemble these multi-level scaling laws and take into account 
critical brain dynamics. 
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Decomposition into different spatial scales showed oscillations in a pre-seizure state at 
all levels. Observation of such oscillations in real world data offers characteristics to be 
useful when testing future models. As we showed here, based on critical transition theory, 
the variance's oscillations along with its scaling law are generically characteristics for Hopf 
bifurcation. These results therefore validate previous seizure models assuming a Hopf bifur- 
cation as a main mechanism as the transition point [721 |83l |50]. While the goal in seizure 
prediction is to predict large events, there is growing consensus about the key role played 
by small events, from precursor oscillations to subclinical seizures [SH El]- Future models 
and predictor systems should encompass those as prediction algorithms unable to account 
for such small oscillations would be ill-adapted and likely provide incorrect seizure forecasts. 

Near seizure onset we observed a time shifted increase in phase-locking; phase-locking 
for higher frequencies (lower scales) tended to precede lower frequencies (higher scales). In 
a recent study, wavelet analysis of spike-wave discharges, a different form seizure activiy, 
revealed changes in the time-frequency dynamics during discharges. While initially a short 
period with the highest frequency value was observed, the frequency later decreased [71 [8] . 
Other studies showed high frequency oscillations specifically at seizure onset [SZl , see [HZ] 
for a comprehensive overview. Together these studies demonstrate dynamic changes in the 
time-frequency domain of seizures with higher dominating frequencies at seizure onset. Our 
observations complement these findings suggesting a frequency-dependent, shifted start of 
synchronization near seizure onset. 



A Data 

Eight patients undergoing surgical treatment for intractable epilepsy participated in the 
study. Patients underwent a craniotomy for subdural placement of electrode grids and strips 
followed by continuous video and electrocorticogram (ECoG) monitoring to localize epilep- 
togenic zones. Solely clinical considerations determined the placement of electrodes and the 
duration of monitoring. All patients provided informed consent. The study protocols were 
approved by the Ethics Committee of the Technical University Dresden. ECoG signals were 
recorded by the clinical EEG system (epas 128, Natus Medical Incorporated) and bandpass 
filtered between 0.53 Hz and 70 Hz. Data were continuously sampled at a frequency of 200 
Hz (patients 1 — 3 and 5 — 8) and 256 Hz (patient 4, [23]). We always indicate the sampling 
point number on the time axis if we use the data. No claims regarding a large-scale statistical 
validity of the data set is made since the total patient sample size is rather small. Although 
this is an important issue [77] we focus here on identifying the dynamical mechanisms and 
new time series analysis techniques in the context of epileptic seizures. 
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